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Our previously developed Constrained-Pairing Mean-Field Theory (CPMFT) is shown to map 
onto an Unrestricted Hartree-Fock (UHF) type method if one imposes a corresponding pair con- 
straint to the correlation problem that forces occupation numbers to occur in pairs adding to 1. In 
this new version, CPMFT has all the advantages of standard independent particle models (orbitals 
and orbital energies, to mention a few), yet unlike UHF, it can dissociate polyatomic molecules to 
the correct ground-state restricted open-shell Hartree-Fock atoms or fragments. 
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I. INTRODUCTION 

In a recent series of papers^— we have devel- 
oped constrained-pairing mean- field theory (CPMFT), a 
method capable of describing static (strong) correlation 
in an accurate and efficient manner. The idea behind 
CPMFT is to make use of the pairing correlations (see 
below) that occur in a quasiparticlc picture to describe 
static correlation in molecular systems. In CPMFT, we 
divide the natural orbitals into core, active, and virtual 
blocks; each core orbital has unit occupation, each vir- 
tual orbital has zero occupation, and the active natural 
orbitals have fractional occupations rij , where < n^ < 1 . 
Static correlation is introduced by allowing electron pairs 
to have fractional occupations within an active space. 

The use of a pairing interaction has many advantages. 
Unlike unrestricted Hartree-Fock (UHF), CPMFT has 
zero spin density everywhere for closed-shell systems. In 
the absence of static correlation, CPMFT reduces to re- 
stricted Hartree-Fock (RHF), while it dissociates poly- 
atomic molecules to restricted open shell Hartree-Fock 
(ROHF) atoms or fragments. Essentially, the dissocia- 
tion limit of CPMFT can be thought of as an ensemble 
solution. By reducing to RHF in the absence of strong 
correlation and ROHF at dissociation, CPMFT cleanly 
separates static from dynamic correlation, as previously 
shown in Ref. y, where the CPMFT P and K density ma- 
trices were used to construct alternative densities to be 
used as inputs into traditional density functionals for the 
dynamical correlation energy. Remarkably, CPMFT ac- 
complishes these feats at a mean field computational cost 
instead of the combinatorial blowup of complete active 
space (CASSCF) or full configuration interaction (FCI). 

While CPMFT is clearly distinct from UHF, it shows 



some unexpected connections. We can take advantage 
of these connections to simplify the formalism, make it 
more efficient, and establish interesting similarities. The 
purpose of this paper is to demonstrate these relations 
and the accompanying reformulation of CPMFT. Ac- 
cordingly, we discuss this connection in Sec. IIIII at some 
length, and show how we can use it to simplify the so- 
lution of the CPMFT equations. Section iPVl shows some 
numerical examples, and we provide conclusions in Sec. 
IVl We include an Appendix that discusses some other 
formal properties of CPMFT. First, however, we provide 
a brief introduction to pairing correlations in Sec. |Hj 



II. PAIRING CORRELATIONS AND THE 
QUASIPARTICLE PICTURE 

Strong correlations in nuclear physics or superconduc- 
tivity are often described as the formation of Cooper 
pairs. The theoretical machinery which docs this is the 
Hartree-Fock-Bogoliubov (HFB) method^. In HFB, we 
write the wave function I'&hfb) as a single determinant 
of quasiparticles created by quasiparticle creation oper- 
ators which are linear combinations of electron creation 
and annihilation operators. The quasiparticle wave func- 
tion thus violates electron number conservation. Because 
the quasiparticle wave function is a single determinant, 
its associated density matrix R is idempotent (R 2 = R) 
and Hermitian (R = TV). We have 



R 



7 

— K* 



K 

1-7 J 



(1) 



Here, 7 is the physical density matrix in the spinorbital 
basis; it is Hermitian but not idempotent. Information 
about pairing correlations is carried by the anomalous 
density matrix k, which is antisymmetric by definition 



because k* 
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We limit our discussion to the 



closed shell case, in which case wc have, for aa, a/3, j3a, 



and /3/3 blocks 



7 
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(2a) 
(2b) 



where P is the closed-shell (spatial orbital) density ma- 
trix and K is the (symmetric, positive semi-definite) 
closed-shell anomalous density matrix. We emphasize 
here that only the a/3 and /3 a blocks of n are non-zero, 
so that we consider only singlet pairing^ We should also 
mention that the notation here differs slightly from that 
used in Refs. [O-y, but does so in an attempt to make 
this manuscript self-contained and as clear as possible. 

Idcmpotcncc of the quasiparticlc density matrix R 
yields two conditions on the electronic density matrix P 
and the anomalous density matrix K: 



PK-KP=0 


(3a) 


P P 2 = K 2 . 


(3b) 


the "odd-electron 


distribution" of 



Physically, K 2 is 
YamaguchifS. the "density of effectively unpaired elec- 
trons" of Starovcrov and Davidson^ and is related to 
Mayer's "free valence index"— once it is written in terms 
of the total density matrix -y aa +7 /3 ^ i = 2P. Essentially, 
K 2 gauges the singlet diradical character of the system 
(or, for larger active spaces, the polyradical character) 
and is a local measure of electron entanglement. 

The HFB energy is given as the expectation value of 
the Hamiltonian with respect to the HFB wave function 



E 



HFB 



= <* 



HFB 



H\$ 



HFB/ 



— 1h P 



(2{ij\kl) - (ij\lk))P lk P 3l 



(4a) 
(4b) 



(ij\kl)KijK kl 



where summation over repeated indices here and 
throughout the manuscript is implied; hij are matrix el- 
ements of the one-electron part of the Hamiltonian and 
(ij\kl) are two-electron integrals in Dirac notation. 

In order to determine the occupation numbers and nat- 
ural orbitals, HFB variationally minimizes £hfb subject 
to the constraint that the density matrix P contains the 
correct number of particles: 



Tr(P) = N. 



(5) 



This condition is enforced by a chemical potential fi in- 
troduced as a Lagrange multiplier. The HFB formulation 
leads to equations similar to Hartrce-Fock, which in the 
particular case of closed-shell systems are 



R cs ?i 



HFB 



u 



HFB 



R r 



0, 



(0) 



where R cs is the closed shell quasiparticlc density matrix 

' K 1 - P ) ^ 



and "Hhfb is the doublc-Hamiltonian (DH) given by 



n 



HFB 



F cs + fiN 
A 



A 

-F cs - liN 



(8) 



Here F cs is the standard closed-shell Fock matrix and A 
is known as the pairing Hamiltonian. These are given by 



F% = hij + (2(ik\jl) - (ik\lj))P kl , 



Aij = (ij\kl)K k i. 



(9a) 
(9b) 



The double-Hamiltonian "Hhfb is just the mean-field of 
the physical Hamiltonian with respect to the quasiparti- 
clc determinant. 

Because the pairing energy of HFB (the term propor- 
tional to K 2 in Eqn. |4b|) is positive when the electron- 
electron interaction is repulsive, the variationally opti- 
mal solution is always K = and therefore A = 0. In 
other words, HFB just returns the regular Hartree-Fock 
solution for Coulombic repulsive systems. In order to 
have HFB solutions with energies lower than Hartrec- 
Fock, one needs a net attractive two-body interaction, 
as in the Bardeen-Cooper-Schriefer picture of supercon- 
ductivity (where it is provided by electron-phonon cou- 
pling) or in nuclear forces. In order to take advan- 
tage of the pairing picture for the conventional repulsive 
electron-electron interaction, and with the aim of describ- 
ing strong correlations, CPMFT simply reverses the sign 
of the pairing energy. We thus have 



E, 



CPMFT 



zn^j ±ij 



(2(ij\kl) - (ij\lk))PikPji 



(ij\kl)KijK k i 



(10) 



The last term plays the role of a correlation energy - 
a correction to the closed shell RHF-like energy expres- 
sion - and will be referred to as such throughout this 
manuscript, but it is certainly not our previous defini- 
tion of static correlation^ -Ecpmft — Ekhf, since P is 
not Prhf- 

In addition to changing the sign of the pairing energy, 
in CPMFT we also restrict non-integer occupations to an 
active space, so that pairing only occurs between quaside- 
generate orbitals. Changing the sign of the pairing term 
changes the sign of A so that the double-Hamiltonian is 



n 



CPMFT — 



■UN 



-A 
-F cs - fiN 



(11) 



Otherwise, DH-CPMFT follows the same procedure as 
in HFB. However, changing the sign of the pairing en- 
ergy and the pairing matrix severs the connection be- 
tween the HFB wave function |$hfb) and the CPMFT 
energy. Note that what we have called simply CPMFT 
in Refs. [[HI is here referred to as DH-CPMFT, whereas 
"CPMFT" here refers to the new formulation to be in- 
troduced below. 

We can, indeed, view the CPMFT energy as the ex- 
pectation value of a model Hamiltonian with respect to 



a particle-number violating determinant: 



H m = K 



CPMFT 



*>, 



,t. 



H = (F$ + h ij )aja j 

2 lJ ~ 



(12a) 
(12b) 



-Ay-aja] - -A^-OiOj. 



This quadratic model Hamiltonian, however, is 7io£ the 
mean-field of the physical Hamiltonian with respect to a 
quasiparticle determinant. As previously noted^ we can 
interpret the CPMFT energy as a hybrid of Hartree-Fock 
and HFB where Hartree-Fock uses 2/7*12 as the electron- 
electron repulsion operator and HFB uses — 1/7*12 • 

Nevertheless, we have a fruitful alternative viewpoint, 
which is to envision the CPMFT energy expression of 
Eqn. (|10p as defining a model two-particle density matrix 
r such that the energy in the spin-orbital basis is 



-Ecpmft — Tr(h7) + Tr(vrcpMFT) 



(13) 



where v is the two-particle part of the Hamiltonian, and 
h is the one-particle part. In terms of spin-orbitals, we 
have 

(Tcpivift)^ = -% 7j - -£ii 7j - 2 i j *• ' 

with lower (upper) indices corresponding to bra (ket) 
indices. The first two terms in this model two-particle 
density matrix correspond to Hartree-Fock whereas the 
last term introduces static correlation via K, which is 
a measure of non-idempotency for P. This last term is 
an important quantity in the cumulant decomposition 
of density matrices^ but in our work appears naturally 
from the idempotency of the quasiparticle density matrix. 
If we use this model two-particle density matrix to de- 
fine expectation values of two-particle operators, then as 
shown in the Appendix, we find the important result that 
CPMFT has no particle number fluctuations. In making 
this choice, we are inevitably working with a density ma- 
trix functional and are effectively doing some form of a 
statistical ensemble theory. Table fl] collects results about 
the UHF two-particle density matrix, the CPMFT model 
two-particle density matrix, and the analogously defined 
HFB model two-particle density matrix. We derive these 
results in the Appendix. 

Note that our model two-particle density matrix 
has appeared before in the literature as the corrected 
Hartree-Fock (CHF) functional^ Our model is, however, 
solved by diagonalization£ More importantly, CPMFT 
restricts the non-integer occupation numbers to an ac- 
tive space only, and in the present work further enforces 
the corresponding pairs constraint, which we will now 
introduce. 



III. CPMFT AND UHF 

Consider the UHF treatment of a system where the 
number of spin-up and spin-down electrons is the same. 



TABLE I: Summary of properties in UHF, HFB, and CPMFT 
for closed-shell systems. We show the correlation energy [i.e 
the difference between the energy from the method and the 
closed-shell piece), the effective polarization, and the particle 
number fluctuations. 



Method 



E c a 



Polarization 






UHF 
HFB 

CPMFT 



-t# M\ M{ 



kl j^i] 
kl RV 



V{- K 13 K kl 



-v 



ij 



1,1 



M = (A-B)/2 

K= |A-B|/2 4Tr(K 2 ) 
K= |A-B|/2 



a v^ = (ij\kl) is a two-electron integral in Dirac notation 



The spin-up and spin-down density matrices 7™ and 
7^ are both idempotent: 



(■y aa ) 2 - 7 aQ = (7 /3/3 ) 2 - 7 /3/3 



0. 



(15) 



The charge density and spin magnetization (or polariza- 
tion) density matrices are 



M = i(7 QQ -7' 9 ' 3 ). 



(16a) 
(16b) 



Traditionally, the UHF energy^ is expressed in terms of 
the ~f aa and 7 /3/3 density matrices: 

£uHF = M7<r+7f/) (17) 

+ \(m)(i? k a +iT)(it +^1) 

where we have put 7°" and 7" in the same basis (say, 
the atomic orbital basis). Although it is almost never 
presented in this way, we can also write the UHF energy 
as a functional of P and M, which yields 

£ UHF [P, M] = £ CS [P] + E C [M], (18a) 

£ CS [P] = 2h ii P ii (18b) 

+ (2(ij\kl) - (iJllk^PikPji 

E C [M] = -(ij\kl)M u M jk . (18c) 

Here, E cs indicates the usual RHF energy expression 
given in terms of the charge density matrix P , while E c 
carries the correlation energy in terms of the spin mag- 
netization density matrix M. An utterly unexpected re- 
sult is that the closed-shell CPMFT energy expression of 
Eqn. (JTUJ) is identical to the UHF energy expression of 
Eqn. (p~8|) . except that the spin density matrix M is re- 
placed by the anomalous density matrix K J^ In cases in 
which UHF predicts static correlation by breaking sym- 
metry (i.e non-zero spin density)^ P is not idempotent. 



Instead, it satisfies 



P P 2 = -(j aa + 7 /3 ' 3 ) - -(~f aa + 7^) 2 (19a) 



= \h aa -l f3IJ ? 
= M 2 . 



(19b) 
(19c) 



This is one consequence of the idempotence of 7 Qa and 
7' 3 ^. The second is 



PM + MP = M. 



(20) 



Note that the condition of Eqn. (fT9|) is the same as the 
CPMFT condition of Eqn. ([3b]). again with M taking 
the role of K. Both the magnetization density matrix M 
and the anomalous density matrix K are Hermitian. 

While CPMFT and UHF thus use the same energy 
expression (one with K and the other with M), K and 
M are not the same even though with the same density 
matrix P, we have K 2 = M 2 . There are also some other 
important differences. Both UHF and CPMFT impose 
an additional condition on these two matrices, which in 
UHF is given in Eqn. $ZU§ while in CPMFT is instead 
given in Eqn. pap . Additionally, K is positive semi- 
definite while M is traceless (and thus has both positive 
and negative eigenvalues). Finally, because in UHF we 
write P as the half-sum of two idempotent matrices, its 
eigenvalues occur in what is known as "corresponding 
pairs" n% and 1 — n ^ 14 ' 15 a terminology that we here 
adopt. 

That UHF has the corresponding pairs property has 
little to do with UHF per se. It originates simply from 
the observation^ that the eigenvalues of a matrix that 
is the half-sum of two idempotent matrices are 0, 1, -^, 
or a corresponding pair (n, 1 — n). Similarly, the eigen- 
values of a matrix written as the half-difference of two 
idempotent matrices are 0, i^, or a corresponding pair 
(— n,n)^' Thus, for example, M has eigenvalues adding 
to in pairs while P has eigenvalues adding to 1 in pairs. 
Quite generally, any non- integer eigenvalues of the charge 
density matrix from a single determinant method will be 
either i or occur in a corresponding pair. Eigenvalues 
of \ could be part of a corresponding pair (for entan- 
gled electrons) or may occur singly for open shells. We 
should be clear that while matrices written as the sum of 
two idempotent matrices exhibit the corresponding pairs 
property, the converse is not necessarily true; a matrix 
whose eigenvalues come in corresponding pairs may or 
may not be the sum of two idempotents. 

Unlike UHF, the eigenvalues of P in DH-CPMFT do 
not occur in corresponding pairs (except when the ac- 
tive space consists of two spatial orbitals). That said, 
the corresponding pairs property has some attractive fea- 
tures for CPMFT. Most important is that it eliminates 
overcorrclation between orbital pairs in different symme- 
tries. This is ubiquitous for example in N2 where the 
variational principle drives occupancy into orbitals at low 



energies and one must introduce multiple chemical po- 
tentials to retain the correct total number of a and 7r 
electrons. A corresponding pair constraint controls this 
unphysical "spilling" and has the inherent attractive fea- 
ture of limiting strong correlations to be an affair between 
orbital pairs. 

Previously, we had introduced the corresponding pairs 
feature within the DH-CPMFT framework using differ- 
ent chemical potentials (Lagrange multipliers) for differ- 
ent irreducible representations of the system. However, 
in the general case where no spatial symmetry is present, 
imposition of this constraint leads to one Lagrange multi- 
plier per orbital pair and a rather complicated nonlinear 
optimization problem. A more satisfactory and much 
simpler approach, however, is to write the CPMFT den- 
sity matrix as 



P = i(A + B) 



(21) 



where A and B are auxiliary density matrices, individ- 
ually idempotent and Hermitian (A 2 = A = A^ and 
similarly for B). As with UHF, the decomposition above 
enforces the corresponding pairs condition automatically, 
and there is no need to enforce this condition via La- 
grange multipliers. Eigenvalues of or 1 in P cor- 
respond to virtual or core orbitals, respectively, while 
paired eigenvalues correspond to active orbitals. Fur- 
ther, by choosing A and B to trace to half the number 
of electrons, we guarantee that P does likewise, and we 
thus have no need of any chemical potential. By making 
this decomposition, in other words, we can avoid the La- 
grange multipliers of the double-Hamiltonian approach 
entirely, and thus simplify the computation. Note that 
once we have converged solutions for A and B (and thus 
P and K), we could, if desired, extract the Lagrange 
multipliers of the DH-CPMFT approach. 

The critical mathematical difference between CPMFT 
as formulated in this manner and UHF is that in UHF, 
we get M from the spin-up and spin-down density ma- 
trices, while in CPMFT, we get K from the total den- 
sity matrix alone (since K satisfies the condition of Eqn. 
(|3b[) . commutes with P, and is positive semi-definite). In 
other words, CPMFT with corresponding pairs defines P 
from A and B as in Eqn. (f2"Tj) , but differs from UHF in 
constructing 

K=v / P^P 7 =iv / (A-B)^I|A-B|. (22) 

from auxiliary density matrices A and B while UHF 
builds P and M from 7"" and 7^, as shown in Eqn. 
(TT51) . Note in the last equation our definition of the abso- 
lute value of a matrix from the square root of the square. 
In practice, to calculate the absolute value of a matrix one 
needs to diagonalize it, flip the sign of the negative eigen- 
values and transform back to the original basis. Both the 
square root and absolute value of a matrix are positive 
definite matrices and both have a convergent polynomial 
series expansion if the matrix is positive definite with 
eigenvalues between and 1, as is the case here. 



To make the comparison between CPMFT and UHF 
more concrete, consider the case where A and B are 2x2 
matrices and let M = |(A — B). Idempotency of A and 
B requires that in the natural orbital basis we have 



(23a) 



B = 



P = 



M 



K 




k = \Jn{l — n )- 



(23b) 



(23c) 



(23d) 



(23e) 
(23f) 



When A and B arc of larger dimension, then in the natu- 
ral orbital basis they are block diagonal with 2x2 blocks 
of the form given above. This is essentially a consequence 
of Eqn. (j2"0")l . which in the natural basis becomes 



(m + n,j)Mi. 



Mi. 



1 



(A, 3 - By), (24) 



the solutions to which are Mi 3 
Because we also have An + Bi 



and rii + rij = 1. 

2 m 5ij , we conclude 
that for i ^ j, we must either have A i3 = By = or 
rii + rij = 1 (in other words, the two eigenvalues form a 
corresponding pair). When the occupation numbers are 
degenerate, the natural orbitals are not uniquely defined 
and we can thus choose them such that A and B still have 
this structure. In the core (virtual) space, A = B = 1 
(A = B = 0). 

Before we continue to the working equations for 
CPMFT in this UHF-likc framework, let us pause to 
make it explicit that CPMFT and UHF are different 
methods. While we have expressed the UHF energy as 
a density matrix functional, we could also write it as an 
expectation value 



E 



UHF 



(<!> 



UHF 



ff|$ 



UHF 



(25) 



with |$uhf) constrained to be a single determinant. This 
is not true of the CPMFT energy expression, and in 
fact there seems to be no wave function associated with 
CPMFT. This may seem somewhat surprising, in light 
of the intimate connection between CPMFT and HFB 
theory, in which there certainly is a wave function, al- 
beit one which violates particle number conservation. As 
we have said, we lose the HFB wave function because we 
have by fiat changed the sign of the pairing energy. Addi- 
tionally, unlike UHF, the spin density is zero everywhere 
for closed shells, even in the presence of static correlation. 
One might wonder whether CPMFT is equivalent to 
projected UHF (PUHF). It is not. If one projects 



the UHF determinant onto a spin eigenfunction, one 
finds that the charge density matrix of the UHF 
determinant and the spin-projected state have the 
same eigenf unctions J^ Spin projection, in other words, 
changes only the occupation numbers of the charge den- 
sity matrix, but not the natural orbitals. The fact 
that the UHF and CPMFT natural orbitals are differ- 
ent should lay to rest any concerns that CPMFT is just 
a projected UHF. 

Another fundamental difference between CPMFT and 
UHF is the onset of the appearance of the solution with 
energy lower than RHF. As shown in our previous paper ja 
the CPMFT solution for a two-level model system ap- 
pears when the RHF orbital energy gap reduces to 

£2 - £i < i<ll|ll> + i<22|22) + (11|22), (26) 

whereas the UHF Coulson-Fischcr instability point is de- 
termined by 



e 2 -si < (12|12) + (11|22). 



(27) 



Because all two-electron integrals in the equations above 
are positive, the CPMFT solution appears inevitably 
when the orbital gap closes and strong correlation is man- 
ifest, such as along a dissociation curve. 



A. Working Equations 

Let us now return to the solution of the CPMFT equa- 
tions in this UHF-like framework. For convenience, we 
repeat the energy expression here: 



E, 



CPMFT 



= E r 



{ij\kl)KijK kl . 



(28) 



We simply minimize the energy with respect to (idempo- 
tent) A and B matrices. The derivatives of E cs in Eqn. 
(|28| with respect to A and B give the usual closed-shell 
Fock matrix obtained from P. That is 



dE c 
dA, 



dE a 

dB, 



ldE c 



2 dPi, 



(29) 



The differences with UHF arise from differentiating 
the last term of the CPMFT energy, which we shall call 
^CPmft_ taking derivatives with respect to A leads to 
an effective potential A, given by 



A, 



0KCPMFT 0KCPMFT Q K 



dA,. 



dK k 



dA,, 



-2A 



dK 



kl 



kl 



dA, 



(30) 
This is essentially the same result that we get from dif- 
ferentiating £™ F of Eqn. OF 



dff c UHF 



dE™dM kl = mF 0M kl 

dM kl d^ a ~ kl 8^ a 



where 



aUHF 
^kl 



(km\nl)M mn — {kl\mn)M„ 



(31) 



(32) 



looks just like A except we replace K with M. In UHF, 
however, we simply have 



dM k i 



1 



ij Z 



SikSji 



(33) 



while in CPMFT the derivative of K with respect to A is 
obtained by differentiating both sides of K 2 = j(A— B) 2 . 
This gives 

^^K ml + K k J^± = \ (M 3l S kl + M hl S 3l ) . (34) 



1 ',/ 



L o 



In the natural orbital basis where K is diagonal with 
eigenvalues Ki, we have 



dK k i _ 1 MjiSkj + M ki Sji 
8A ll ~ 2 K k + Ki 



(35) 



Thus, in the natural orbital basis the effective potential 
A is 



A,;, = - 



AuMji A kj M k 



Ki + Ki K k + Kj 



Since 



dK 
dA 



dK 
9B' 



(36) 



(37) 



the equations we ultimately solve are [F A ,A] = and 
[F B , B] = 0, where F A and F B are effective Fock matrices 
given by 



F A = F c 
F B = F c 



A, 

A. 



(38a) 
(38b) 



At first glance, the right-hand-side of Eqn. [36] might 
appear to be divergent unless all K t are non-zero. How- 
ever, since forcing A,-j = actually gives the condition 
Kij = 0, we simply set A,-j = for the inactive- inactive 
(core and virtual) block where K must be zero (because 
the occupation numbers are or 1). Therefore, in Eqn. 
1361 such divergent terms due to inactive orbitals are sim- 
ply removed from the sum. 



IV. RESULTS 

We have implemented this version of CPMFT in the 
Gaussian suite of programs^ Each calculation requires 
the specification of the number A^Act of active natural 
orbitals. Due to the corresponding pairs constraint, the 
number of active electrons is always equal to iVAct ~ in 
other words, we always work at half-filling. In order to 
obtain an appropriate initial guess for A and B, we mix 
the coefficients of the AAct orbitals closest to the Fermi 
level, just as one would do to break spatial symmetries 
in UHF. The natural orbital pairs closest to the Fermi 
energy correspond to those whose occupations are closest 
to half and half. 



In single bond systems where we normally choose the 
active space to be two electrons in two orbitals, the cor- 
responding pair constraint is automatically satisfied, and 
no difference is observed between the results using the 
present approach and those using our previous double- 
Hamiltonian approach (that is, diagonalization of the 
doublc-Hamiltonian constructed from F and A). How- 
ever, in DH-CPMFT, one must adjust the chemical po- 
tential /x at every iteration of the SCF procedure to con- 
trol the number of electrons in the active space. Be- 
cause we must adjust the chemical potential, we must 
diagonalize the double Hamiltonian of Eqn. (fTTj) sev- 
eral times in each SCF cycle, until the resulting density 
matrix has the proper trace. In contrast, the current 
approach requires no chemical potential, since we have 
Tr(P) = 1/2 Tr(A + B). Because both A and B trace to 
the correct number of electrons, so too does P. This is 
a significant operational advantage of the present imple- 
mentation. 

For systems with larger active spaces, the present ap- 
proach differs from DH-CPMFT, although as mentioned 
above, we can impose the corresponding pairs constraint 
in DH-CPMFT in some special cases by including differ- 
ent chemical potentials for different irreducible represen- 
tations. We illustrate this with the case of N2. Table [TT1 
shows the total energy of N 2 at 2.0 A. We use the cc- 
pVTZ basis set and choose six active orbitals and six ac- 
tive electrons. The current scheme gives a slightly higher 
energy than does DH-CPMFT with only one chemical 
potential, as one would expect since we have imposed an 
additional constraint on the system. Also as one would 
expect, it gives the same results as does DH-CPMFT 
with the corresponding pairs constraint enforced by ad- 
ditional Lagrange multipliers. However, removing the 
chemical potentials results in considerable computational 
savings. In Fig. [TJ we show the N2 dissociation curves 
from CPMFT in the doublc-Hamiltonian approach and 
in the corresponding pairs framework. In this case, the 
corresponding pairs constraint has only a minor effect on 
the energy. 

We have also performed a CPMFT calculation of the 
C2 molecule with the 6-31G basis set. Near equilib- 
rium, C2 has significant static correlation due to near- 
degeneracy between the RHF occupied a\ s and unoccu- 
pied G2 Pz orbitals. As the molecule is stretched, how- 
ever, the ir x , TT y , 71-*, and ir* orbitals become degener- 
ate, while the c^-o^ interaction becomes weak. We 
have therefore chosen our active space to be six electrons 
in six orbitals for this system. In Fig. [5] we show the 
total energy of C2 as a function of bond length. The 
CASSCF energy includes all static correlation that re- 
sults from these orbital interactions (plus some dynam- 
ical correlation). Without the corresponding pairs con- 
straint, DH-CPMFT strongly overcorrelates nearly ev- 
erywhere. Adding the corresponding pairs constraint sig- 
nificantly reduces this overcorrclation. Near equilibrium, 
it gives results between UHF and CASSCF. Unfortu- 
nately, it still overcorrelates as the molecule dissociates. 



TABLE II: CPMFT energies of N2 at R = 2.0 A. Also included are the number of diagonalization steps required, Adiag, and 
the number of SCF cycles required for convergence. 



Scheme 



Energy (a.u.) N dia 



DH-CPMFT(6,6) a 
DH-CPMFT(6,6) 6 
CPMFT(6) 



-108.79901762 
-108.79715442 
-108.79715442 



"Single chemical potential 

'Corresponding pairs enforced by multiple chemical potentials 



118 

121 
12 



SCF cycles 



32 

34 

12 



-108.7 
-108.75 

-108.8 
-108.85 

-108.9 

-108.95 

-109 

-109.05 

-109.1 
-109.15 



0.5 




RHF 

UHF 

CASSCF(6,6) 

DH-CPMFT(6,6) 

CPMFT(6) 



2.5 




RHF 

UHF 

CPMFT(4) 

CPMFT(6) 

CASSCF(6,6) 

DH-CPMFT(6,6) 



1.5 



Rn-n (A) 



2 2.5 

Rc-c (A) 



3.5 



FIG. 1: Potential energy curves of N2 calculated with the FIG. 2: Potential energy curves of C2 calculated with the 
cc-pVTZ basis set. 6-3 1G basis set. 



This is due to electron "spilling" between a% s and <72p z or- 
bitals. As R — > 00, only the 7r orbitals should be strongly 
correlated; including these a orbitals in the active space 
at large internuclear separation allows them to correlate 
and lower the energy unphysically. If we remove two or- 
bitals from the active space, we produce the curve marked 
CPMFT(4). This goes to the correct dissociation limit, 
but undercorrelates at equilibrium where the active space 
should be larger. The correct solution for this molecule 
involves introducing renormalized one-body potentials in 
CPMFT(6) that eliminate the spilling at dissociation, 2 
an approach that we will discuss in a forthcoming article. 
While going to the right dissociation limit is important, 
it is perhaps less critical than getting the correct behav- 
ior near equilibrium. Note that CPMFT (4) dissociates 
correctly to two ROHF carbon atoms, while UHF in- 
stead dissociates to two spin-contaminated UHF carbon 
atoms and CASSCF(6,6) has some dynamical correlation 
at dissociation. 

Finally, we stress the differences between UHF and 
CPMFT by analyzing the dissociation of the CO2 
molecule. The ground state of CO2 near equilibrium is 
a closed-shell singlet with no expected static correlation. 
Indeed both UHF and CPMFT reduce to the RHF so- 
lution near R e . However, when the molecule is symmet- 
rically stretched and the two oxygen atoms arc simul- 
taneously separated from the carbon atom, the correct 
dissociation limit corresponds to all three atoms in their 



triplet ground state. This situation cannot be handled 
by UHF. In CO2 near R e , there are six electrons associ- 
ated with bond formation, three with spin-up and three 
with spin-down. At dissociation, UHF might assign two 
spin-up electrons to one oxygen atom and two spin-down 
electrons on the other, which puts both oxygen atoms 
in their triplet ground state. However, with only one 
electron of each spin remaining, the best UHF can do 
is to assign a singlet state to the carbon atom, which is 
clearly incorrect and not the lowest energy state. In sim- 
ple words, UHF runs out of broken symmetry degrees of 
freedom (has only two) to model the dissociation of CO2 
(Fig. [3]) and misses the correct dissociation limit by ~ 
20 milliHartrecs. The bumps in the dissociation curves 
correspond to crossings of different solutions to the re- 
spective SCF equations and we have plotted the lowest 
energy state at each R. Because spin states are treated in 
CPMFT through an "ensemble" representation, one that 
yields zero spin magnetization density everywhere, the 
CPMFT solution for this dissociation has two half spins 
up and two half spins down on each of the three atoms, 
leading to the correct energy corresponding to the sum of 
ROHF atomic energies. Note that CPMFT(6) in Fig. [3J 
contains a one-body potential arising from an asymptotic 
constraint as explained in our previous publication.— We 
defer detailed discussion of the renormalization schemes 
used in CO2 and applicable to C2 within the current 
UHF-like context to a forthcoming publication. 




RHF 

UHF 

CASSCF(6,6) 

CPMFT(6) 

ROHF Fragments. 

UHF Fragments 



2.5 



one mapping between A and B on the one hand and P 
and K on the other. At dissociation, for example, solu- 
tions where A and B orbitals are localized and delocal- 
ized (roughly corresponding to UHF and RHF orbitals) 
are degenerate. The existence of additional degenerate 
solutions in CPMFT (compared to UHF) can lead to 
convergence difficulties as the active space becomes large. 
Efficient ways of dealing with the additional degrees of 
freedom provided by the auxiliary A and B matrices are 
currently under investigation. 
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FIG. 3: Potential energy curves for the double dissociation of 
CO2 calculated with the 3-21G basis set. 



V. CONCLUSIONS 

We have developed a novel scheme for performing 
CPMFT calculations with occupation numbers occur- 
ring in corresponding pairs. In doing so, we eliminate 
all chemical potentials, and the effective Fock matrices 
F A and F B that are to be diagonalized are of half the 
dimension of the double Hamiltonian matrix in the pre- 
vious DH-CPMFT scheme. Thus, the computational ef- 
fort in our present implementation is greatly reduced 
over the previous formulation of CPMFT. The corre- 
sponding pairs constraint reduces the overcorrelation of 
C2 near equilibrium, and has important consequences 
for the dissociation of heteronuclear systems. While the 
corresponding pair constraint could also be imposed in 
the DH-CPMFT framework by addition of one Lagrange 
multiplier per electron pair, the current approach im- 
poses this constraint in a simpler black-box manner. 

We have shown that this version of CPMFT is closely 
related to UHF theory. Unlike UHF, however, CPMFT 
incorporates static correlation by a different mechanism. 
The physical density matrix 7 has identical spin-up and 
spin-down blocks, whereas the auxiliary A and B density 
matrices, in general, break symmetry. CPMFT can cor- 
rectly dissociate polyatomic molecules into ROHF atoms 
or fragments, whereas UHF has problems with multiple 
entangled electrons at multiple centers, as shown for CO2 
above. In the present formulation, CPMFT becomes a 
density matrix functional that can be solved by diagonal- 
ization of effective Fock matrices providing orbitals and 
orbital energies. We wish to emphasize one more time 
that as we have demonstrated, a quasiparticle picture of 
strong correlations with the sign of the pairing interac- 
tion reversed yields an energy expression reminiscent of 
UHF. 

Finally, we should note that in CPMFT different auxil- 
iary A and B density matrices can lead to solutions with 
degenerate energies. The key quantities determining the 
energy in the model are P and K and there is a many-to- 
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Appendix A: Properties of the CPMFT Model 
Two-Particle Density Matrix 

The CPMFT model two-particle density matrix is 



(T 



CPMFT J 



\kl 



1 



iH 



1 



H 



(Al) 



where i, j, k, and I are spin-orbitals and 7 and k are the 
density matrix and anomalous density matrix in the spin- 
orbital basis (i.e, they are of dimension 2 A x 2 A, where 
A is the size of the atomic orbital basis). In general, 7 is 
Hermitian and k is antisymmetric. When everything is 
real (which we take for simplicity; this does not affect our 
conclusions), the idempotent HFB quasiparticle density 
matrix is 




R 



Idcmpotcncy tells us that 

■y k — n-y = 0, 

2 2 

7 — K = 7. 

We recall that for closed shells £ 




(A2) 



(A3a) 
(A3b) 



(A4a) 

(A4b) 

(A4c) 

(A4d) 



We can define an analogous model two-particle density 
matrix for HFB, for which all the conditions on k, 7, K, 
and P are the same, but where 



(Phfb) 



1 



.k-.l 



1 



.l.Jc 



1 



kl 



lili + n^jK 



K-l 



(A5) 



Finally, the UHF two-particle density matrix is 



(1 UHF )ij - -jli lj - -7t7j 


(A6) 


where 7 is idempotent. We have 




(~1 aa \ /P + M \ 

7 ~ ^ j^J ~ [ p - mJ ' 


(A7a) 


P = P 2 + M 2 , 


(A7b) 


M = PM + MP. 


(A7c) 



1. Partial Trace of the Two-Particle Density 
Matrix 



An important condition on the two-particle density 
matrix is that it traces to the one-particle density ma- 
trix. That is, we must have 



v >y 



(A8) 



We remind the reader that repeated indices arc to be 
summed. 

The partial trace condition is satisfied by the UHF two- 
matrix and the CPMFT model two-matrix, but not by 
the HFB model two-matrix: 



2 Wi - 7i7J =F KijK a ) 



^[TV 7 <-( 7 2 )<±(* 2 ) 



-[TV 7 <-( 7 + * 2 )<±(k 2 ) 



^7< - \ [(« a )$ T (* 2 )'] 



(A9a) 
(A9b) 
(A9c) 
(A9d) 



2. Particle Number Fluctuations 

In order to work out particle number fluctuations, we 
need the expectation values of TV and TV 2 , with TV' the 
number operator, given as 



J pq lJ, p u, q- 



(All) 



We have already noted that the expectation value of TV 
is just Tr(7). The expectation value of TV 2 requires the 
two-particle density matrix: 



(TV 2 ) = 8 pq S rs {a\,a q ala s ) 

= 8 pq 5 rs {-{a],a\.a q a s ) + 8 qr {a p a s )) 
= 6 pq S rs (2T^,+S qr %) 



— O-ppr 



ll- 



(A12a) 
(A12b) 
(Al2c) 
(A12d) 



If the two-particle density matrix obeys the partial 
trace condition, the particle number fluctuations are au- 
tomatically zero. This is thus true of UHF and of 
CPMFT. However, HFB has particle number fluctua- 
tions: 



(TV 2 >HFB = (TV-lbj-2(K 2 ^+ 7 j 



TV 2 



2 Tt(k) 



implying that 

a% = (TV 2 ) - (TV) 2 



-2Tr(K 2 ). 



(A13) 



(A14) 



Note that this is positive, as it should be, since ~k 2 = 
7 — 7 and occupation numbers are between and 1, in- 
clusive. In the closed-shell case, we have a 2 N = 4Tr(K 2 ). 



3. Spin Contamination 

Evaluating spin contamination is more complicated 
than evaluating particle number fluctuations, not least 
because we need an expression for (S 2 ) for a general two- 
particle density matrix T. We begin by noting that 



S 2 = S 2 



■si 

SS+, 



(Al5a) 
(A15b) 



Here, the top (bottom) sign in ± and =F corresponds to 
CPMFT (HFB), and we have used antisymmetry of k. 
Explicitly, we have 



(T 



CPMFT 



(r 



TV-l 

TV- 1 



"%• 



HFBj 



""/ 



(« 9 )i. 



(AlOa) 
(AlOb) 



Note that by TV we mean the trace of the one-particle 
density matrix 7, which should be the number of particles 
in the system. 



where S± is the spin raising/lowering operator. We are 
interested here in the closed-shell case (i.e. N a = N/j 
with a block diagonal 7). 

In the closed-shell case, the contribution to (S 2 ) from 
the first term is zero. We must evaluate the contribution 
from the next piece using our model two-particle density 
matrix. We have 



4 2 = X>(*) 2 +I>(^(i) 



(Ai6) 



'f-j 



X- 



The first (second) term is a one-particle (two-particle) 
operator. Note that we could also write 



t = 2^2s z (i)s z (j) 



(A17) 



l>3 



Thus, we end up with 

(Si) 



HFB 



<S 2 )cpmft - Tr(K 2 ) 
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(A24a) 
(A24b) 



which explains the factor of 2 that might otherwise ap- 
pear to be missing below. 

Evaluating the contribution to (S 2 ) from X z is 
straightforward, and we get just 



The contribution to (S 2 ) from S-S+ must also be eval- 
uated using the model two-particle density matrix. Ex- 
panding this operator in terms of contributions from in- 
dividual electrons, we have 



{X z ) = ±(N a +N /3 ) = ±Tr(P). (A18) 



§-§+ = J23-(f)S+(i) + J2$-W+U)- ( A25 ) 



The nonzero matrix elements of Y, are 



x 



(^)i± = (i a J a \Y z \k a l a ) = -6lS{, (Al9a) 

(Y*)i a Z = (iJ(i\£\k a h) = —Sidf, (A19b) 

(Yz)%t = (Wa|£|Ma) = -\ 5 l 5 l ( A19C ) 

(Y*Z% = (i?Jp\Y z \hh) = \s\8l. (A19d) 

Here, we are working in an orthornomal basis set. 

The relevant components of the CPMFT and HFB 
two-particle density matrices are 



r 



\(fi J&- -fa-fa) > ( A20a ) 

\{^l^^j^ kah ), (A20b) 

K^-ytTKi^ 1 "), (A20c) 



,k l a _ 1 kff l a kpl 



■ W/3 



\(fi'if.-j&'&) 



y\kplp J- ' ^Kp^lp 



'lp '3p 'IP <3P 



where the top (bottom) sign corresponds to CPMFT 
(HFB). 

Contracting the density matrices with the matrix ele- 
ments, we get 

(A21) 
where we have used antisymmetry of k. Working in our 
closed-shell case, this reduces to 



(f 2 ) = -iTr(P 2 T K 2 ). 



(A22) 



In total, then, we find that (S' 2 ) in CPMFT and HFB 
is given by 

(S 2 ) = ^Tr(P - P 2 ± K 2 ) (A23a) 



-Tr(K 2 ± K 2 ) 



(A23b) 



The first term is the one-particle operator X, and the 
second is the two-particle operator Y. 

Since X does nothing to down-spin electrons but an- 
nihilates up-spin electrons, we clearly have 



(X) =N = Tr(P). 



(A26) 



To take the expectation value of Y, it proves useful to 
symmetrize it so that it acts the same on the two elec- 
trons. Since operators acting on different electrons com- 
mute, we have 



f = £M*)Mj) 



(A27a) 



'¥./ 



i£(M*)Mj) + M^-0-)) (A27b) 



It-j 



53 (S_(t)s+(j) + *+(*)*-(?)) 



l>] 



The only nonzero matrix elements of Y are 



Y 



ipla 



Y, 



ia3P 

kal a 



{ipja\Y\kJ p ) =S l k 6l 
{i a jp\Y\k p l a ) =S i k 6 J l . 



(A27c) 



(A28a) 
(A28b) 



The relevant spin components of the CPMFT and HFB 
model two-particle density matrix are 



r 






-p,fc Q Lfi 



Ic, k P 

lie Ijp 

Jp ,k a 



I ( h k a k a \ 



(A29a) 
(A29b) 



where again CPMFT (HFB) corresponds to the top (bot- 
tom) sign. 

Contracting the two-particle density matrix with the 
matrix elements gives us 



(Y) = -Tr(7 aQ ,7^g =F MV)' 



(A30) 



In the closed-shell case, using the results in Eqn. (|A4|) . 
this becomes 



(Y) = -Tr(P 2 =F K 2 



(A31) 
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Then the expectation value of S- S+ is given by 

(£_£+) =Tr(P-P 2 TK 2 ) 
= Tr(K 2 TK 2 ). 

We therefore have 

(S-S + )kfb = 2Tr(K 2 ), 

(S-S+ )cPMFT = 0. 



(A32a) 
(A32b) 

(A33a) 
(A33b) 



For UHF in cases in which there is strong correlation, 
we have the familiar formula 



(S 2 ) = s(s + 1)+Nf3- Tr( 7QQ 7w ). (A35) 



For the closed-shell case, using the results in Eqn. (|A7jl . 
we have 



Combining Eqns. (|A24|) and (|A33|) gives us the total 
spin contamination in HFB and in CPMFT: 



(5 2 )hfb = 2Tr(K 2 ), 
<£ 2 )cpmft = Tr(K 2 ), 



(A34a) 
(A34b) 



(S 2 ) = Tr[P - (P + M)(P - M)] 

= Tr(P-P 2 +M 2 ) 
= 2Tr(M 2 ). 



(A36) 
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